Simplified Onsager theory for isotropic-nematic phase equihbria of length 

polydisperse hard rods 
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Polydispersity is believed to have important effects on the formation of liquid crystal phases in 
suspensions of rod-like particles. To understand such effects, we analyse the phase behaviour of 
thin hard rods with length polydispersity. Our treatment is based on a simplified Onsager theory, 
obtained by truncating the series expansion of the angular dependence of the excluded volume. We 
describe the model and give the full phase equilibrium equations; these are then solved numerically 
using the moment free energy method which reduces the problem from one with an infinite number 
of conserved densities to one with a finite number of effective densities that are moments of the 
full density distribution. The method yields exactly the onset of nematic ordering. Beyond this, 
results are approximate but we show that they can be made essentially arbitrarily precise by adding 
adaptively chosen extra moments, while still avoiding the numerical complications of a direct solution 
of the full phase equilibrium conditions. 

We investigate in detail the phase behaviour of systems with three different length distributions; 
a (unimodal) Schulz distribution, a bidisperse distribution and a bimodal mixture of two Schulz 
distributions which interpolates between these two cases. A three-phase isotropic-nematic-nematic 
coexistence region is shown to exist for the bimodal and bidisperse length distributions if the ratio 
of long and short rod lengths is sufficiently large, but not for the unimodal one. We systematically 
explore the topology of the phase diagram as a function of the width of the length distribution and 
of the rod length ratio in the bidisperse and bimodal cases. 



I. INTRODUCTION 

Suspensions of rod-like particles can undergo an orien- 
tational disorder-order phase transition from a disordered 
isotropic (I) phase, in which rods point with equal prob- 
ability in every direction, to an orientationally ordered 
nematic (N) phase in which particle orientations cluster 
around a preferred direction. This behajuiaur has been 
obsecyed experimentally in both chemicaU □ and biolog- 
icalDuij systems. 

Two main theoretical approaches have been used to e: 
plain the I-N phase transition. The Maier-Saupe theor 
neglects density variations and is therefore appropriate 
for so-called thermotropic materials, where the transi- 
tion is driven by variations in temperature. It focuses on 
attractive, long-range interactions between particles; in 
the original theory these attractions were thought of as 
arising from Van der Waals forces but in fact the Maier- 
Saupe approach can be considered a lowest order approx- 
imation fat. general orientation-dependent attractive in- 
teractionsEj. The Onsager theoryEJ, on the other hand, 
takes into account only hard core repulsions between the 
particles on contact; temperature then becomes unimpor- 
tant and so the theory is appropriate for lyotropic ma- 
terials, where the phase ordering is driven by changes in 
density. The I-N phase transition arises from the compe- 
tition between the orientational entropy (corresponding 
to the tendency of rods to stay orientationally disordered) 
and the packing entropy (which is higher for aligned rods 
due to the excluded volume interaction). The key simpli- 
fication is that, in the "Onsager limit" of infinitely thin 
rods, the virial expansion for the free energy truncates 
after the second-order contribution. The state of the 



system is described by the density of rods and the dis- 
tribution P{9) of the angles 9 which the rods make with 
the nematic axis, i.e., the preferred orientation in the ne- 
matic phase. A formal minimization of the free energy 
with respect to this function leads to a self-consistency 
equation for P{9). Solving this gives, in principle, the 
free energy as a function of density, from which the I- 
N phase coexistence region can then be obtained by a 
standard double-tangent construction. Onsager gave an 
approximate solution for the I-N phase transition using 
a variational ansatz for P{9). With two jdifferent appropt 
imation techniques Lekkerkerker et al.E3 and Lashertj 
obtained more accurate resultp^the numericaJly exact so- 
lution has also been obtaineal3'Lj. IsiharaEj, basing his 
calculations essentially on the Onsager theory, had shown 
earlier that the order-disorder phase transition will occur 
at similar densities for particles of different shapes, as 
long as they are strongly asymmetric (e.g. plates, cylin- 
ders, spheroids). 

We mention briefly an alternative to Oiisager theoty 
that was developed by Floryll3, Di MarzicO and AlbenEa, 
based on a lattice model. The interaction is again of ex- 
cluded volume type, but particles are now embedded in 
a three-dimensional cubic lattice whose sites are occu- 
pied either by segments of the rod-like particle or by sol- 
vent. Although qualitetively this model gives the same 
behaviour as the Onsager theory, the I-N phase transi- 
tion is pieedicted at a density about twice that fopid by 
Onsagertata. With a similar lattice model FloryE3 pre- 
dicted an analogous order-disorder phase transition for 
semiflexihlfi macromolecules. On the experiinental side, 
RobinsonEll, HermansEj and Nakajima et alE found re- 
sults in good agreement with Flory's theory for hard rod- 
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like molecules. 

Returning now to Onsager theory, a direct compari- 
son between the theoretical predictions and experimen- 
tal data is difficult due to the complexity of experimen- 
tal systems. For example, there arp forrections due to 
non-hard interactions in real systemstJO which affect the 
densities of coexisting isotropic and nematic nhascs and 
other characteristics of the phase transitionE3cj. These 
effects can be partially accounted for by adding an aiiuac.- 
tive potential to the excluded volume interactionEZlc3, 
which leads to a temperature dependence of e.g. the ori- 
entational order parameter. 

More importantly for our purposes, polydispersity - 
the presence of an (effectively) continuous spread of rod 
lengths and/or diameters - has been|-sjbe.wn to haAte-iim- 
portant effects on the phase diagramcBE^I. OnsageJlJ al- 
ready gave the extension of his approach to the case of 
mixtures of rods with different lengths or different diam- 
eters, but the obvious complication of having many (and 
in the fully polydisperse limit, infinitely many) distin- 
guishable particle species prevented a direct numerical 
solution, or the systematic use of trial functions as in the 
monodisperse case. We focus in this paper on length poly- 
dispersity alone, assuming that rod diameters are still 
monodisperse. Such length polydispersity can lead to 
a strong broadening of the coG^astence region, which has 
been observed experimentallyE3c3 and also predicted the- 
oretically, atjleaat within an expansion for narrow length 
distributionsEilS. Coexistence of several nematic phases 
(N-N), occasionally with a third isotropic phase (I-N-N), 
has also been seen experimentallycJ and been predicted 
by theory for bi- and tridisperse systems (comprisi|rig.r|Ods 
of two and three different lengths, respiec1p:cly)ll3'LJTJ. 

In a series of papers by Flory et al.EJL^, the length- 
polydisperse case was solved within the lattice model for 
a bidisperse system as well as for cases with many dif- 
ferent rod lengths (which in the lattice model are always 
multiples of a given unit length). The results obtained 
are qualitatiiiely similar to the ones predicted by the On- 
sager theorylij, with a three-phase I-N-N region and an 
N-N region at high density for the bidisperse system, and 
with strong fractionation and broadening of the coexis- 
tence regioitior length distributions of e.g. PoissonEj and 
exponential type. Moscicki et al. applied the same the- 
ory to a system with a Gaussian length distributionO, 
observing strong fractionation and a pronounced depen- 
dence of the size of the coexistence region on the width 
of the distribution. 

The theoretical studies referred to above which con- 
cern themselves directly with the off-lattice situation are 
restricted to either mixtures of rods of two or three differ- 
ent lengths or to continuous length distributions of small 
width. The question of the effects of full length polydis- 
persity within the Onsager theory of thin hard rods there- 
fore remains open, and it is this problem that we address 
in the present paper. It is clear that the introduction of 
full polydispersity significantly complicates the analysis 
of phase behaviour. Each rod length now has its own 



conserved density associated with it, and a naive appli- 
cation of the Gibbs phase rule suggests that one could 
have arbitrarily many coexisting phases. This problem 
is compounded by the presence of the varying rod orien- 
tations, so that the full density distribution over lengths 
and orientations has both conserved and non-conserved 
parts. A direct attack on the Onsager model with full 
length polydispersity is therefore extremely difficult. 

As a first step towards understanding the phase be- 
haviour of polydisperse hard rods, we have previously 
analysed the Zwanzig model, a discretization of the On- 
sager theory which allows rcpds, to point only in one of 
three orthogonal directionsLJc3. We found that this ex- 
hibited some of the features that are typical of polydis- 
perse systems (broadening of the coexistence region, frac- 
tionation of rod lengths across different phases). It was 
not, however, able to predict the observed N-N and I-N-N 
phase separations. This could be traced back to the re- 
striction on rod orientations: a separation of a single ne- 
matic phase into two nematics containing predominantly 
shorter and longer rods involves a loss of entropy of mix- 
ing and a gain in orientational entropy; with only three 
different rod orientations available, the possible gain in 
orientational entropy is too limited for such a phase sep- 
aration to be favourable. 

In the present work we analyse a model which can take 
the complexities of full length polydispersity into account 
and still predict the more complex phase coexistences in- 
volving several nematic phases that were missing in the 
Zwanzig model. Our "V2 Onsager model" is an approxi- 
mation of the full Onsager theory obtained by truncating 
the expansion of the angular dependence of the excluded 
volume of two rods after the lowest nontrivial term, which 
involves the Legendre polynomial 1^2 ■ The model is thus 
no longer exact, even in the limit of thin rods, but nev- 
ertheless useful as a solvable approximation to the full 
Onsager theory which can, for example, provide some in- 
formation about the possible topologies of the phase di- 
agram. The errors introduced by the truncation can be 
directly assessed for the bidisperse case, where accurate 
numerical solutions of the full Onsager theory existE£lE3, 
and we find below that there the V2 Onsager model re- 
produces most of the trends in the evolution of the phase 
diagram with the rod length distribution. The truncation 
errors could in principle be reduced by extending the an- 
gular expansion to increasingly high orders, giving a se- 
ries of models whose behaviour will approach that of the 
full Onsager theory in the limit; the V2 Onsager model is 
the first member of this series. Its free energy also bears a 
formal similarity to that of the polydisperse Maier-Saupe 
model, and by imposing the constraint of constant den- 
sity could yield qualitative information about the phase 
behaviour of polydisperse thermotropics. A final benefit 
is that the V2 Onsager model is, in contrast to the full 
Onsager theory, truncatable: its excess free energy only 
depends on a finite number (two) of moments of the den- 
sity distribution. Phase equilibria can therefore be fepnd 
efficiently using the moment free energy methodE3E3 (or 
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moment method for short) , with its extension to systems 
with non-conserved densities introduced in the context 
of the polydisperse Zwanzig modeled. 

This paper is structured as follows. In Sec. |l| wc de- 
scribe our truncated "7^2 Onsager model" and its link to 
the full Onsager theory, and give the exact phase coex- 
istence equations for the truncated model. In Sec. Ill, 



we briefly review the moment method and then apply it 
to the V2 Onsager model. The numerical method used 
for calculating phase equilibria is described in some de- 
tail, including our procedure for ensuring that the mo- 
ment method produces numerically accurate predictions 
for phase equilibria even beyond the point where it is ex- 
act by construction (which turns out to be the onset of 
nematic phase ordering). Our main results are detailed in 
Sec. where we give phase diagrams for three different 
types of rod length distributions: Schulz distributions, 
which are unimodal, bidisperse distributions containing 
only two different rod lengths, and finally mixtures of 
Schulz distributions that interpolate between these two 
extremes. One of our guiding questions will be under 
which conditions on the rod length distribution the phase 
diagram shows the more "exotic" features of nematic re- 
entrance and of I-N-N and N-N coexistence regions. We 
will find that, while these features are absent for the uni- 
modal Schulz distributions with even the largest widths 
that we consider, the bidisperse and bimodal systems give 
a much richer phase behaviour. We summarize in Sec. 
and discuss directions for future work. App. ^ contains a 
description of a simple approximation that we found use- 
ful in elucidating certain trends in the evolution of the 
phase diagram, particularly in the limit of large length 
ratios between long and short rods for the bimodal and 
bidisperse distributions. 



II. MODEL DETAILS 

We consider a system of hard cylinders of length L and 
diameter Z?, capped with hemispheres at both ends. We 
assume that these rod-like particles are polydisperse in 
their length L, but all have identical diameter D. To 
be able to take the Onsager limit of long thin rods, we 
introduce a reference length Lq and the normalized rod 
lengths I — L/Lq; the Onsager limit is then reached for 
D/Lq — > at fixed values of I . 

The state of a single phase of a system is described 
by the density distribution p{l,il), defined so that 
p{l , ^) dl dft / 4:TT gives the density of rods with (normal- 
ized) lengths in an interval dl around and orientations 
O in a solid angle dfi. The distribution p{l, fl) is the nat- 
ural extension for a polydisperse system of the usual ori- 
entational dijatribution P{Q) used for monodisperse sys- 
tems of rodafJ. The rod orientation can be parame- 
terized in terms of the angle 6 with the nematic axis and 
an azimuthal angle if] due to the cylindrical symmetry 
of the nematic phase p(/, fl) is independent of ip. Using 



dil = d cos 9 dip the density distribution as a function of 
rod length, obtained by integrating over orientations, is 
therefore 



p(i) = 1- I dn p{i,a) = I de p{i,i 



Here and below we use the shorthand 



de F{0) 



d cos 61 F{0) 



(1) 



(2) 



for the angular integral of an arbitrary function F{9). 
The orientationalj-distribution of rods can be factored 
out from p(l, 9) asE^ 



p{i,e)^p{i)Pi{e) 



(3) 



where Pi [0) represents the probability of finding a rod of 
given length I in orientation O — (9, ip) and is normalized 
to 1: 



d9 Pii9)^l 



(4) 



In the isotropic phase, one has Pi{9) = 1 and p{l,9) = 
p{l). 



A. Polydisperse Onsager theory 

As mentioned above, the present model is a trunca- 
tion of the full Onsager theory (see e.g. Refs. 32 4^ for 
reviews), which we now describe. In the Onsager limit 
D/Lq — > 0, the second order virial approximation be- 
comes exact and gives for the excess free energy in units 
of fceT 



/ = Jdi dl' ^ ^ B2{iLo,n,no,n')p{i,9)p{i',e') 

(5) 

Here B2 is the second virial coefficient which for hard 
rods is half the excluded volume and given by 



B2 = DLl Z/'|sin7| 



(6) 



up to neghgible terms of order D/Lq, with 7 the angle 
between the two rods. Using the cylindrical symmetry of 
the nematic phase one can perform the integration over 
Lp and ip' in Eq. (H). If we also make the densities / and 
p{l, 9) non-dimensional by multiplying by a unit volume 
{tt/A)DLq we obtain 



/ = - t dldl' d9 d9' W K(9, 9') p(l, 9)p(l', 9') (7) 

TT 



where the angular part, 
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K( 



2^27"'^"^" 



A;2„7'2n(cos6')p2n(cos 9') 



n=0 



(8) 



can be expcessed as a bilinear expansion in Legendre 
polynomialalJ 7'2n(cos0) with fco ~ 7r/4, fc2 — — 57r/32 
etc. and 7'o(cos6') = 1, V2{cose) = i(3cos2 6'- 1). The 
excess free energy is therefore an infinite sum of moments 
of the density distribution p(l,0), with weight functions 
/p2n(cos6»): 



/ 



ra=0 



k2n dl dO I V2nicos 9) p{l, 9) 



(9) 



B. The 7^2 Onsager model 

From Eq. (|) we see that the full Onsager theory does 
not give a truncatable excess free energyo, i.e. one that 
is just a function of a finite number of moments pt of 
the density distribution p{l,9). To obtain a more man- 
ageable-Lheory, to which the efficient moment free energy 
methodoCJ can be applied, we now truncate the excess 
free energy after the second Legendre polynomial. We 
will thus consider the "7^2 Onsager model" defined by 
the excess free energy 



/ 



Cl 2 



C2 2 

yP2 



where 



and 



pi= dl d9 lp{l,9) 



P2^ dl d9 lV2icos 9) p{l, 9) 



(10) 



(11) 



(12) 



are moments of the density distribution with weight func- 
tions 



'Wi{l,0) = I 

W2{l,e) = lV2icos9) 



(13) 
(14) 



and the numerical constants are ci = (8/7r)fco = 2, 
C2 = — (8/7r)fc2 = 5/4. In our units pi is the rescaled 
volume fraction (j) of rods, pi = {Lo/D)(l>. We will denote 
below by po the zeroth moment of the density distribu- 
tion, po = J dlp{l), with weight function wo{l,9) = 1; 
this is the number density of rods. The ratio pi/po = (l) 
then gives the rod length averaged over the normalized 
length distribution P{1) = p{l)/po, and P2/P1 < 1 ex- 
presses the degree of nematic ordering. 

It is clear that the truncated model defined above will 
give approximate results compared to the predictions 



of the full Onsager theory. However, previous workM 
for the monodisperse case has investigated the conver- 
gence of the densities of coexisting isotropic and nematic 
phases when the expansion (|^) is truncated at higher and 
higher orders. Already for the truncation after 1^2 the re- 
sults were qualitatively correct, and this encourages us to 
study polydispersity effects within the truncated model 
defined above. 

A different rationale for studying the 7^2 Onsager 
model could come from the fact that its excess free en- 
ergy (|lO|) is very similar to that of the polydisperse Maier- 
Saupe model for thermotropics (see e.g. Refs. ^2|). There 
the interaction is not of the hard core type, but if one 
assumes that the interactions of rods of lengths I and I' 
again scales as W the final expression for the excess free 
energy is identical to Eq. ( [To| ) , with ci = and C2 = C2 (T) 
a function of temperature. The actual predictions of the 
two theories would nevertheless not necessarily be close, 
since within Maier-Saupe theory the overall rod density 
Pa is assumed to be the same in all phases; this implies, 
for example, that a finite region of I-N coexistence can- 
not appear for monodisperse rods (see e.g. Ref. We 
therefore leave a study of polydisperse Maier-Saupe the- 
ory for future work. 

To complete the specification of our model, we need 
to add to the free energy the ideal part, i.e. the free en- 
ergy of an ideal mixture of polydisperse particles (see e.g. 
Refs. 11,0,51 



/id= / d9 dl p{l, 9) [In p{l, 9)- 1] 



(15) 



Using the decomposition (^) and the relation (||), the 
total free energy (density) of our model is therefore: 



I dl p{l) [\np{l) - 1] + I dldO p{l)Pi{9) \nPi{9) + / 

(16) 



with / given by ([lOD. For a given density distribution 
p{l), Pi{9) is obtained by minimization of the free en- 
ergy with respect of Pi{6), subject to the constraint (|^). 
Introducing corresponding Lagrange multipliers k{1) one 
has the condition 



= 



SPii9) 



/+ dl d0 k{1)Pi{9) 



(17) 



= p[l) [\nPi{0) + 1] + lp{l) [cipi - C2P2V2] + nil) (18) 
Solving for Pi{0) gives (here and in the following 7^2 = 

P2(C0S(?)) 



Pi{0) 



exp(Zc2/02'P2) 

/ d0 exp(Zc2P2'P2) 



(19) 



The orientational distributions are therefore determined 
only by p2, which from ( p^ ) obeys the self-consistency 
equation 
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P2= [ dl dO p{l)- ^ 

J J d9 exp(/c2P2'P2) 



(20) 



If there are several solutions P2 for a given p{l), the one 
with the smaller free energy (Jig) is the physical one. 



C. Phase coexistence equations 

At this point we need the expressions for the chemical 
potential fi{l) and the (osmotic) pressure 11 in order to 
derive the phase equilibrium conditions. The chemical 
potential is obtained by functional differentiation of the 
free energy ( p^ ) with respect to p{l); the variations of 
the Pi{0) with p{l) need not be considered since Pi{9) is 
chosen to minimize /. After a little simplification one 
thus finds 

^(0 = \np{l) -\njd9 exp{-lcipi + IC2P2V2) (21) 

The pressure follows from the Gibbs-Duhcm relation as 

n = -f + j dl p{l)p{l) = po + yP? - ^pI (22) 

Imposing equality of chemical potentials in a set of P 
coexisting phases, labelled by a = 1 ... P, we can express 
the density distributions as 



^^-^^(0 = R{1) J dO exp(-kip[") + lc2pt'V2) (23) 

where R{1) is a function of I common to all phases. Using 
Eq. (^^, the full density distributions over rod lengths 
and orientations follow as 

p^''\l,e)^R{l)eM~lcip[^^ +lc2p["^V2) (24) 



The function R{1) can be found from Eq. ( P3|) and the re- 
quirement of particle conservation: if p^"-* (Z) is the overall 
or "parent" density distribution and the fraction 

of the system volume occupied by phase a, then 



This leads to 
R{1) 



j2v^'^^pi^Hi)-P^"Hi) 



p(o)(Z) 



vi-) J d9 eM-lcip[^^ + lc2pt'V2) 



(25) 



(26) 



so that, from Eqs. ( |l9| , p3| ), the density distributions in 
the coexisting phases are 



p(")(0exp(-?cipr+^C2Pr^2) 



(a) 



/ dO eM-lcJi^ + IC2P'2^V2 



(27) 



In principle we thus have the conditions for coexistence 
between P phases: In each phase, p2 obeys the self- 
consistency condition (pO|) and pi a similar equation 



obtained by inserting ( p7| ) into (|l]); for the P phase 
volume fractions w^"^ we have the equality of pressure 
in all phases (P — 1 equations) plus the normalization 
w^"^ — 1. However, finding a starting point from 
which a numerical solution of this strongly coupled sys- 
tem of nonlinear equations converges successfully is very 
difficult. We use instead the moment free energy method, 
which gives exact results for some features of the phase 
diagram and otherwise allows us to approach the solution 
of the full phase equilibrium conditions with controllable 
accuracy. 



III. THE MOMENT METHOD 

We only outline the construction of the moment free 
energy here and refer to Ref. ^ for details of its prop- 
erties, and Ref. 43 for the extension to non-conserved 



degrees of freedom such as the rod orientations in our 
case. Terms which are linear in the conserved densities 
p{l) = J dO p{l,0) can be added to the free energy ( [l^ ) 
without affecting the phase equilibria, since they merely 
add constants to the chemical potentials p{l). We are 
therefore free to replace / by 



/ = de dl pH, 



r(0 



/ 



(28) 



with r(l) an arbitrary function of I. Guided by the intu- 
ition that it is the moment densities pi and p2 appearing 
in the free energy that drive phase separation, we now 
allow violations of the lever rule (^) as long as they do 
not affect these moments. All other details of the density 
distribution p(l, 6) are then found by minimizing the free 
energy (Ea) for the given values of pi and p2 as defined 
in Eqs. (OJT2). Inserting a Lagrange multiplier for 
each moment pi {i ~ 1,2), the minimum value of the free 
energy is then 



mom ~ ^ ^ ^iPi 



PO + f 



(29) 



where po = J dl p{l) is the number density of rods as 
defined previously. The density distribution at which this 
minimum free energy is obtained is 



p{l, 9) = r{l) exp ^ \wi{l, 



(30) 



where the Wi(Z, 9) are the weight functions (|l^,|l^ defin- 
ing the moments; the are determined implicitly by the 
requirement that the density distribution ( ^0|) gives the 
correct values for the pi. 

Eq. ( p9| ) defines the moment free energy] rather than 
being a functional of the density distribution p(l), i.e. of 
an infinite number of conserved densities, it is a func- 
tion of the densities pi and p2- To find phase equilibria 
using the moment free energy, one proceeds as follows. 
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The "moment chemical potentialsji-jare defined as pi = 
dfnwm/dpi and can be written aai3 pi = Xi + df /dpi, 
giving in our case 



Hi = Ai 



cipi 



= A2 - C2P2 



(31) 



(32) 



Since pi is conserved while p2, which contains the rod 
orientations, is not, Hi has to be identical in a set of 
coexisting phases while 112 has to vanish in all the phases. 
The pressure calculated from /mom, 

n = -/mom + Pi/ii + P2m = Po + yP? - y P2 (33) 

also has to be identical in all phases. Finally, the lever 
rule has to be satisfied for the conserved moment pi, i.e. 

(0) (a) {"■) 

It is now easy to see that the moment free energy gives 
the correct phase equilibrium conditions, apart from the 
violations of the lever rule that it allows. In fact, equality 
of pi and vanishing p2 imply from Eq. ( |30| ) that the 
coexisting phases calculated from the moment free energy 
can be written as 

p^'^) (1,6)^ r{l) exp (/ii? - Icipi""^ " lc2pi"^'P2) (34) 

This is of exactly the form ( p^ ) that we derived earlier 
from the full phase equilibrium conditions and the mini- 
mization of the free energy w.r.t. the rod angle distribu- 
tions. The requirement of equal pressures in all phases is 
also captured exactly by the moment free energy, since 
the expression ( ^3| ) actually coincides with the one de- 
rived from the original free energy of the V2 Onsager 
model, Eq. (H). 

We have seen so far that any phase equilibria cal- 
culated from the moment free energy obey the exact 
conditions of equality of chemical potentials and pres- 
sure. The lever rule will be satisfied only for pi but 
not, in general, for the whole of p{l). However, this fi- 
nal requirement will clearly also be fulfilled if all phases 
in an exactly calculated phase split of the 7^2 Onsager 
model are of the form (|3^). We can guarantee that 
this is the case at least until the onset of nematic or- 
der (the so-called "cloud point" of the isotropic phase), 
by choosing the so far unspecified function r{l) to be 
equal to the parent density distribution p'^°)(Z). This 
is the choice we make from now on. It works because 
at the isotropic cloud point the nematic phase still oc- 
cupies a negligible fraction of the system volume, so 
that only the isotropic phase contributes to the denom- 
inator of Eq. (p7|). Explicitly one gets from Eq. (27) 



with obvious labels for the phases, p^{l,0) = p*-°''(0 and 
p^{l,9) = p(o)(Oexp[/ci(pl -p^) + lc2P^V2]; as claimed, 

(0)/ 



these are both of the form ( pOD with r{l) = p^''{l). 

Beyond the onset of nematic phase coexistence, phase 
equilibria constructed from the moment free energy will 



in general not solve the full phase equilibrium conditions, 
because the lever rule cannot be satisfied with density dis- 
tributions from the family (pO[). Violations of the lever 
rule can, however, be minimized by enlarging the fam- 
ily ( |30| ) to include the actual density distributions oc- 
curring in the coexisting phases. This can be done by 
retaining extra moments beyond those appearing in the 
excess free energy. Obviously, in the limit of retaining 
infinitely many extra moments, the moment free energy 
would eventually recover the full free energy of the model, 
but any computational efficiency would then be lost. The 
idea is therefore to add only a few extra moments to 
reduce the approximation caused by using the moment 
method, keeping their number sufficiently low to avoid 
making the computation too slow. 

In order to enlarge the family ( ^0| ) we could in princi- 
ple use moments defined by any weight function Wi{l,9). 
Since the corresponding moments pi do not appear in the 
excess free energy, their moment chemical potentials are 
Pi — Xi+df /dpi — Xi. The values of the extra Lagrange 
multipliers must therefore be the same in all phases, and 
from Eq. ( |30| ) we can write the density distributions pre- 
dicted by the moment method as 



pl'^) (?,(?) =p(o)(Z) exp 



X exp 



E 



XiWi{l, 9) 



(35) 



Comparing with Eq. (p7|), one sees that the second ex- 
ponential, in which the extra weight functions appear, 
ought to provide a good approximation to the denomi- 
nator of Eq. ( p^ ) , which is not present in Eq. (|3^) . Since 
that denominator is a function of / only, the same needs 
to be true for the extra weight functions. We will always 
retain the zeroth moment, i.e. the total number density, 
with weight function wq — 1. This ensures that the "di- 
lution line" of density distributions p{l) = const, x p'^^\l) 
is included in the family (|30|), a useful feature when we 
want to scan across the phase diagram in the direction 
of increasing or decreasing density. 

In order to choose the other extra weight functions, 
we apply the recently developed "adaptive method"L3. 
The intuition behind this is as follows. A phase split 
calculated by the moment method for a given point in 
the phase diagram will give a total density distribution 
Ptot(0 = w^°-'p''°HO which in general differs from the 
parent p'^^\l) (though agreeing in the value of the con- 
served moment pi). A convenient measure of this viola- 
tion of the lever rule is the "log-ratio" lnptot(0/P^°HO- 
Since from Eq. (|3^) extra weight functions appear expo- 
nentially in the density distributions, including an extra 
weight function W'i{l) similar in shape to the log-ratio 
should drive yOtot(0 closer to p^^\l)] it extends the fam- 
ily ( ^0[ ) in the "right direction" in order to contain the 
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exact distribution. The method therefore proceeds by it- 
eratively including extra weight functions that are fitted 
to the current log-ratio; a proliferation of extra moments 
is avoided, by combining the extra weight functions appro- 
priatelyH. In this way one never requires more than two 
extra moments in the calculation, but these are contin- 
ually adapted to approach the solution of the full phase 
equilibrium equations. The quality of the aporoximation 
is measured by the average square log-ratioc3 



(^ + 1) 



z+l 



6^ dl In 



Ptot(0 

p(o)(0 



and we could easily reach values of 5 < 10^^, imply- 
ing that our results are essentially identical to those that 
would result from a direct solution of the full phase equi- 
librium conditions for the 7^2 Onsager model. We verified 
this explicitly at some points in the phase diagram: with 
the solution obtained from our procedure as a starting 
point, a Newton- Raphson solver for the full phase equi- 
librium conditions converged and terminated at an es- 
sentially indistinguishable solution. 

To fit the extra weight functions to the log-ratio in 
the above procedure, we initially followed the method 
of Ref. representing the weight functions as linear 
combinations of some fixed set of basis functions. For 
some of the more numerically difficult situations, such as 
bimodal distributions with strongly different weights in 
the two peaks, it turned out to be more efficient instead to 
represent the extra weight functions as spline (piecewise 
cubic polynomial) fits to the log-ratio. 



IV. PHASE DIAGRAMS 

In this section we will present our phase diagrams for 
the length polydisperse 7^2 Onsager model. We consider 
first a Schulz distribution of lengths, which has frequently 
been used to model simple unimodal distributions (see 
e.g. Ref. We will find strong differences compared 
to the results for the bi- and tridisperse (full) Onsager 
model. By contrasting with the predictions of the P2 
Onsager model for bidisperse systems, we show that these 
differences are not due to our truncation of the model but 
rather due to qualitative differences in phase behaviour 
for discrete and continuous length distributions. Finally, 
we study a bimodal mixture of two Schulz distributions. 
This allows us to analyse in detail how the topology of 
the phase diagram changes as we interpolate between the 
two extremes of bidisperse and unimodal distributions. 



A. Unimodal length distribution 

We study a parent phase with density distribution 

P^^HO = Po"''-^'*''' (0 where the normalized length distri- 
bution is of the Schulz form: 



r(z-f 1) 



■Pexp[-(z + l)/] 



(36) 



This gives an average rod length of one for the parent, 
a convention we will follow throughout since a different 
average parent rod length could always be absorbed into 
the reference length Lq. The parameter z controls the 
width of the distribution, more precisely its normalized 
standard deviation cr, which from now on we will refer to 
simply as "polydispersity" : 



{I? 



1 



I 



(37) 



The limit z ^ 00 therefore corresponds to the monodis- 
perse case, where the Schulz distribution degenerates to 
a (5-peak. Decreasing z will give increasing polydisper- 
sity; we limit ourselves to z > and thus cr < 1 since 
otherwise the Schulz distribution ( p6| ) has a power-law 
divergence for l^Q.lTL Fig. I we show, for each value of 




FIG. 1. Rod number densities of the isotropic and ne- 
matic cloud phases (solid) and of the corresponding shadows 
(dashed) , plotted against the polydispersity cr on the vertical 
axis. The isotropic cloud curve has the lower density of the 
two cloud curves and joins the isotropic shadow at cr = 0, 
where the system becomes monodisperse. The nematic cloud 
and shadow curves similarly join in the limit. 



a, the isotropic cloud and nematic shadow - giving the 
rod number densities in the coexisting isotropic and ne- 
matic phases at the point where the nematic first appears 
- and the nematic cloud and isotropic shadow, which re- 
late to the point where the isotropic phase disappears 
as density is increased. As anticipated, the coexistence 
region, which is delimited by the two cloud curves, broad- 
ens strongly with increasing polydispersity. The other in- 
teresting feature, which was not seen in our earlier study 
of polydispersity effects in the Zwanzig modeCl, is the 
crossing between the cloud curves and the corresponding 
shadow curves; for values of a above the crossing points 
the number density of the isotropic phase is larger than 
the one of the coexisting nematic phase. If we switch to 
a different representation of the phase diagram by plot- 
ting the (scaled) rod volume fractions pi rather than the 
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FIG. 2. The scaled rod volume fractions of the isotropic 
and nematic cloud phases (solid) and the relative shadows 
(dashed), plotted against the polydispersity a. 



number densities po in the various phases (Fig. ^ the 
crossing between cloud and shadow curves disappears. So 
the rod volume fraction in the nematic is always larger 
than in the isotropic phase, even though the isotropic 
can have the larger number density. Since pi/po is the 
average rod length in each phase, this is clear evidence 
of a strong fractionation effect, with the longer rods be- 
ing found predominantly in the nematic phase. We show 
this explicitly in Fig. |^. Bearing in mind that the cloud 
phases have the same length distributions as the par- 
ent and therefore an average rod length of unity, one sees 
that indeed the average rod length in the nematic shadow 
is much larger than in its corresponding isotropic cloud, 
and that of the nematic cloud is much larger than that 
of the corresponding isotropic shadow. 
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FIG. 3. Plot of the average rod lengths in the isotropic 
(solid) and nematic (dashed) shadow phases. The average 
rod length in the cloud phases, being identical to that of the 
parent, is always one. Thus the nematic phase always has a 
larger average length than its coexisting isotropic phase. Frac- 
tionation disappears when the system becomes monodisperse, 
o" — > 0, as expected. 



B. Bidisperse rod lengths 

From Figures |l|, || and || we saw that there is no I-N-N 
three-phase region and associated N-N demixing within 
the present model, at least for a Schulz distribution of 
lengths. On the othsp-ihand, previous work on bi- and 
tridisperse systemsE3c3 shows that such a three-phase 
separation is possible within the full Onsager theory. The 
absence of a three-phase region observed above could 
then be due to either the difference in the length dis- 
tributions investigated (continuous and unimodal versus 
discrete with two or three different lengths), or to our 
truncation of the Onsager theory that gave the 7^2 On- 
sager model. To ascertain which of these applies, we now 
study the phase diagram of a bidisperse system of rods 
within the V2 Onsager model. The normalized parent 
length distribution is then 



= [(1 - x)S{l - h) + x5{l - h)] 



(38) 



where li < I2 and x is the (number) fraction of longer 
rods. We again fix the average length to unity, so that li 
and I2 can be written in terms of their ratio q — I2/I1 as 



1 



1 — X + xq 
The polydispersity is given by 

2 x{l-x){q 



I2 



qh 



(1 — a; + xq)'^ 



(39) 



(40) 



and a becomes zero for x = and x — I, where we 
recover a monodisperse system; its maximum for given q 
is reached at a: = l/(g -I- 1). 
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FIG. 4. Phase diagram of a bidisperse system for rod length 
ratio q — 2.5. We only show the limits of the the I-N coex- 
istence region, which are given by the isotropic and nematic 
cloud curves, and omit the shadow curves. All densities shown 
are therefore those of the parent phase. 

In Fig. ^ we show the phase diagranfl for rod length 
ratio q = 2.5. The density of the parent at the tran- 
sition from single-phase to two-phase regions is plotted 
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against the number fraction of longer rods on the verti- 
cal axis; effectively, we are only showing the cloud curves 
and omitting the shadow curves. This "phase boundaries 
only" representation will be useful below when three- 
phase equilibria appear (which are represented on the 
cloud and shadow curves only in terms of the exceptional 
points where direct transitions from one to three phases 
occur). Comparing with the Schulz distribution case in 
Fig. |l], where the solid lines correspond to those of Fig. ^, 
we see that again the coexistence region broadens as we 
move from either of the monodisperse limits x — and 
X — 1 towards larger polydispersity. 

For comparison with previous results on bidisperse sys- 
tems obtained from the full Onsager theorj€3'ca, it will 
also be useful to show our phase diagrams in the represen- 
tation employed in Refs. t 
4)1 = {Li/D)(l)i and 02 
Pq. Here (/>i and (/)2 are the volume fractions of short and 
long rods, and Li = LqIi is the unnormalized length of 
the short rods. With A^i and N2 the number of short and 
long rods, we have 
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which uses the variables 
{Li/ D)(j)2 instead of our x and 



D ' 



D 4 



— -zrnH 



V 



N 



(41) 



The term in brackets is just our dimensionless rod num- 
ber density po, so that 



ii = (1 - x)IIpo 
62 = xhhpo 



(42) 
(43) 



This gives us the relation between the two representa- 
tions of the phase diagram. As an aside, we note that 
our convention for choosing the normalized rod lengths 
li — Li/Lq and I2 — L2/ Lq differs from that of Ref. |3^; 
we vary li and I2 according to Eq. ( ^9| ) to maintain a 
constant average length in the parent, while in Ref. ^ 
constant values of li and I2 are used. However, these as- 
signments can alternatively be thought of as arising from 
a different choice of the reference length Lq for identi- 
cal values of Li and ^2- Since Lq does not enter the 
definitions ( ^ , p3[ ), the values of the 4>i and (j)2 remain 
unaffected. We show the (^1, 02-representation for the 
case q = 2.5 in Fig. ||. 

By including Legendre polynomials up to Vu, and de- 
termining the orientational diS|t«bution functions Pi{9) 
numerically, Lekkerkerker et alx^ obtained good approx- 
imations to the phase diagrams for the full Onsager the- 
ory of bidisperse rods. They found that at q — b there 
is a re-entrant ncmatic phase, while at (7 = |2-|thcrc is no 
such re-entrance. Buining and LekkerkerkcrnSi also found 
no re-entrance at q = 2.5; the latter finding agrees with 
our above results. Using a Gaussian variational approx- 
imation for the PiiO), Odijk and LekkerkerkeiO showed 
that re-entrance is genericajlly to be expected for large q. 
Vroege and LekkerkerkeiH later found, using the same 
approach, that in addition to the re-entrance one obtains 
three-phase I-N-N and two-phase N-N equilibria for q 




FIG. 5. Phase diagram of a bidisperse system with q — 2.5 
in terms of the rescaled volume fractions 01 and 02 of long 
and short rods, respectively. This representation is useful to 
compare with previous results for bidisperse systems from the 
full Onsager theory; it also shows directly at any point of the 
phase diagram the composition of the parent phase. 



above « 3.17. Their approximation predicted that the 
boundaries of the N-N coexistence region should be in- 
dependent of the rod number density po, corresponding 
to horizontal lines in a {po,x) representation of the phase 
diagram or to lines thmugh the origin in a {4>i,4>2) plot. 
Van Roij and MuldertJ later showed that these state- 
ments are exactly true (only) in the limit of large po\ this 
means in particular that the N-N coexistence region|46 
not closed off by an N-N critical point. Birshtein Gl_al.lla, 
using a trial function of the type used by Onsagerllil, had 
earlier come to the opposite conclusion. However, they 
emphasized (as did Abe and FloryLj) the difficulties of 
solving the phase equilibrium conditions at large density, 
which reduced the reliability of their results. 

To compare the above findings from treatments of the 
full Onsager theory with the P2 Onsager model, we now 
investigate how the phase diagram changes with increas- 
ing q. Figs. H, 1^, H and ^ show the phase diagrams for 
rod length ratios q — 3.5, q — 7 and q = 12 in both rep- 
resentations {po,x and 01,02)- At q — 3.5 a re-entrance 
of the nematic phase has almost appeared; at q = 7 this 
is fully developed. At g = 12 (Fig. ^), finally, we have 
in addition a three-phase I-N-N region bordered by a re- 
gion (N-N) of coexistence of two nematicsE2l. The N-N 
region is closed off by an N-N critical point, where the 
moments pi and p2 in the two nematics become identi- 
cal. These features persist at higher q; for q = 30, for 
example, the phase diagram (not shown) has the same 
structure as Fig. |^ but with both the re-entrance and 
the three-phase region being located at smaller x and 
extending to higher densities. For even larger q the nu- 
merical evaluation of the phase diagrams becomes more 
troublesome because the interesting features of the phase 
diagram shift towards smaller and smaller x where even 
with the moment method the phase equilibrium condi- 
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FIG. 6. Phase diagram of a bidisperse system with q = 3.5, 
represented in terms of the parent density po and the long rod 
number fraction x (main plot) and in terms of the rescaled 
volume fractions 0i and (j>2 (inset). The large "bump" of the 
nematic cloud point curve (towards large densities in the main 
figure) narrowly avoids a re-entrance of the nematic phase. 




PO 

FIG. 8. Phase diagram of a bidisperse system with q = 12, 
in terms of the parent density pa and the long rod number 
fraction x. As shown in detail in the inset, a three-phase I-N-N 
region appears, bordered by a region of N-N coexistence; at 
larger values of x, the nematic phase is re-entrant. 



tions cannot be reliably solved numerically. However, a 
simple approximate treatment (see App. ^ can be used 
to make progress and understand how the phase diagram 
scales for large q. 




po 

FIG. 7. Phase diagram of a bidisperse system with q — 7 
in both the {po,x) and the (02, ^ii) representations. The 
re-entrance of the nematic phase is now fully developed, but 
there is not yet an N-N or three-phase (I-N-N) region. 



In summary of our results for the bidisperse case, an 
increase in q produces similar qualitative changes in the 
structure of the phase diagram as for the full Onsager 
theory. The quantitative details of course differ, with 
the onset of both the nematic re-entrance and the three- 
phase region shifted to higher q. The only qualitative 
difference caused by the fact that we are using the trun- 
cated p2 Onsager model is that, in contrast to the full 
Onsager theorycJ, the N-N region does not extend to 
large densities but is instead closed by a critical point, 
whose properties we discuss further in App. H. This dis- 
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FIG. 9. Phase diagram of a bidisperse system with q — 12 
in terms of 01 and 02, the scaled volume fractions of shorter 
and longer rods at q = 12. 



tinction between the P2 Onsager model and the full On- 
sager theory can be made plausible by considering the 
effect of our truncation on the excluded volume term, 
i.e., the excess free energy, in the limit of large densi- 
ties pq where nematic phases will be stropghL ordered. 
In the full Onsager theory, it can be shownoEj that the 
excess free energy density is simply f = 2po for large po 
so that the excluded volume per particle is f/po = 2. 
This simple scaling with po arises because the typical 
rod angles are ~ l/po for large po; the angular part 
K{9, 0') of the excluded volume interaction, as given by 
Eq. (||) , scales linearly with the typical values of 9 and 9' 
and is thus also proportional to 1/po- The high density 
phase behaviour thus arises solely from a competition be- 
tween the two components of the ideal part of the free 
energy (p^, the entropy of mixing and the orientational 
entropy; N-N demixing occurs when the resulting gain 
in orientational entropy dominates the loss of entropy 
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of mixing. In our V2 truncation, on the other hand, K 
is approximated as K{9,9') ^ ci — C2V2{coa9)V2{cos9') 
which becomes ~ a + b[9^ + (6*')^] for small angles, with 
constants a and b. For large po, when the typical values 
of 9 and 9' again become small, this approaches a con- 
stant rather than decrease to zero as in the full Onsager 
theory, giving a quadratic density scaling of the excess 
free energy. The same conclusion can be reached from 
the moment description: for 9 we have p2 = pi and 
the excess free energy density becomes \ {ci — C2)p\- The 
quadratic scaling of the excess free energy with density, 
as opposed to the linear scaling in the full Onsager the- 
ory, will act as a driving force against phase separation. 
This is in line with our findings above regarding the ab- 
sence of N-N demixing at high density. One can in fact 
show formally (see App. ^ that N-N coexistence is not 
possible in the 1^2 Onsager model in the limit of large p^. 

The above argument can in fact be extended to all 
possible truncations of the full Onsager theory, since the 
inclusion of any finite number of higher order Legen- 
dre polynomials will still give a behaviour K{9,9') ~ 
a + b[9^ + {9')'^] of the orientational part of the excluded 
volume for small angles; the correct linear scaling with 9 
and 9' of the full theory is obtained only when all terms 
in the infinite series are retained. Equivalently, in the 
moment description one would have a finite number of 
moments analogous to p2, but these would all saturate 
to values proportional to pi at high density, giving again 
the quadratic scaling f ^ p\. We conclude, therefore, 
that there cannot be any N-N region at sufficiently high 
density for any truncation of the full Onsager theory; tlua 
is consistent with van Roij and Mulder's explanatiorO 
for Birshtein's result of an N-N critical point in his ap- 
proximate treatment of bidisperse Onsager theory. Of 
course, for any finite density we must eventually recover 
the correct phase behaviour as more and more higher or- 
der Legendre polynomials are included to give an increas- 
ingly good representation of K{9^ 9') for nonzero angles. 
One would expect, therefore, that successively higher or- 
der truncations of Onsager theory would lead to N-N 
coexistence regions extending to higher and higher den- 
sities, and eventually diverging as the order of truncation 
is taken to infinity. 

Before passing to a different length distribution, let 
us note that in the volume fraction representation (e.g. 
Fig. ^) of the phase diagrams of our bidisperse system 
all tielines connecting coexisting phases are necessarily 
straight. This implies, in particular, that the boundaries 
of the three phase region, which are formed by two-phase 
tielines, are straight. The same argument no longer ap- 
plies in the scenario that we study next: when more than 
two different rod lengths are present, the two-dimensional 
phase diagrams are cuts through a higher-dimensional 
phase diagrams, and generic tielines no longer lie within 
this cut plane. 



C. Bimodal length distribution 

Our results for the bidisperse case suggest that the V2 
Onsager model correctly reproduces most of the quali- 
tative features of the phase diagram of the full Onsager 
theory, except for N-N coexistence at high densities. This 
suggests similar qualitative agreement also in the case 
of unimodal length distributions. We therefore conclude 
that the differences between unimodal and bidisperse dis- 
tributions that we found, in particular the absence of re- 
entrant features and N-N and I-N-N phase equilibria in 
the unimodal case, are physical and not due to our trun- 
cation of the model. To understand in more detail how 
these features arise in dependence on the shape of the 
length distribution, we now turn to a family of length dis- 
tributions that lets us interpolate between the two limits 
of bidisperse and unimodal distributions. The normal- 
ized parent length distributions we consider are mixtures 
of two Schulz distributions, of the form 

(0 = (1 - x)l^^S (l/h) + xl^^S (l/h) (44) 

Here S{1) denotes the Schulz distribution (|3^) and h and 
I2 are again given by Eq. ( |39| ) to constrain the average 
rod length to unity. Let us call cts the normalized width 
(standard deviation) of the Schulz distribution, to distin- 
guish it from the normalized overall standard deviation 
a of the parent distribution; the distribution (^) is then 
characterized by the three parameters q, x and erg. For 
CTs ^ 0, we recover the bidisperse case; on the other 
hand, as as increases for given x and g, the distribution 
eventually becomes unimodal as shown in Fig. For 
given q and x, the value of as at which this happens is 
determined by the condition that there must be a value of 
I for which P^°'>{1) has vanishing first and second deriva- 
tive; the resulting set of equations can easily be solved 
numerically for <ts. 

Our aim is now to understand how the occurrence of 
the more exotic features of the phase diagram, such as an 
I-N-N three-phase region or a re-entrant nematic phase, 
depends on the rod length distribution; one might, for 
example, suspect a dependence on bimodality (the pres- 
ence of two separate peaks in the distribution) or on the 
overall polydispersity. 

We will proceed in two separate ways. First, we fix a 
value of q (=12) for which in the bidisperse case (us — 0) 
we have a three-phase region and re-entrance in the ne- 
matic phase. We then increase erg ^-nd monitor the cor- 
responding variation of the (po, x) phase diagram. When 
as becomes sufficiently large, we have a unimodal length 
distribution and thus expect to see behaviour more akin 
to that for the Schulz length distribution. For fixed q, 
the three-phase region should thus disappear when erg 
reaches a sufficiently large value. The process will essen- 
tially be the opposite of what we have done so far: in the 
previous sections we began from a unimodal phase dia- 
gram and moved to the bidisperse case with larger and 
larger rod length ratio q. Now we start from the largest 
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FIG. 10. Examples of our bimodal parent length distribu- 
tion. In the limit where the width of each of the two peaks 
tends to zero, the bidisperse case is recovered. As the width 
is increased, on the other hand, the distribution eventually 
becomes unimodal; in the inset, we show the case where the 
minimum and the second maximum have just merged into a 
turning point. 



value of q that we used for a bidisperse distribution, and 
widen the two peaks of the distribution until we return 
to a unimodal (though not a simple Schulz) distribution. 

The second procedure will be to fix values of x and 
q for which we have coexistence of three phases in the 
bidisperse limit as 0, and to start increasing the poly- 
dispersity. This lets us check whether the existence of a 
three-phase region is linked directly to the bimodality of 
the parent, by comparing the value of as for which the re- 
gion of three-phase equilibrium disappears to that where 
the parent length distribution changes from bimodal to 
unimodal. 

Adopting the first procedure, we have found the phase 
diagrams for a small normalized width as = 0.1 of the 
two peaks in the bimodal length distribution. For q — 3.5 
the result is a phase diagram that is practically identical 
to Fig. ^; the system is still effectively bidisperse. More 
interestingly, the same conclusion applies when we in- 
troduce a small peak width for q = 12: with as = 0.1, 
we again get phase diagrams that are essentially indis- 
tinguishale from the bidisperse case shown in Figs. ^ 
and ^. In particular, neither the three-phase region nor 
the re-entrant nematic phase are significantly affected by 
such a small degree of polydispersity of the two peaks. 
Larger values of ag, on the other hand, do have a sig- 
nificant effect. If the width of the peaks is increased 
to as = 0.6, for example (Fig. 0), the phase diagram 
changes noticeably. In particular, the three-phase region 
has now disappeared, or at least shrunk to values below 
X = 10~^ where our numerical calculations become un- 
reliable, and the nematic re-entrance has become much 
less pronounced. This is in line with the intuition that for 
increasing as behaviour similar to that for the unimodal 
Schulz distribution should be recovered. The results so 
far confirm that, in order for a three-phase region to ex- 
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FIG. 11. Phase diagram for a bimodal length distribution 
with rod length ratio q — 12 and normalized width of the 
two peaks of as — 0.6. Compared to the bidisperse case with 
the same q (Fig. ^ , the three-phase (I-N-N) region as well as 
the associated N-N region have disappeared. The re-entrance 
of the nematic phase is still present, but is now visibly less 
pronounced than in the bidisperse case or for small as, e.g. 
as = 0.1. Towards a; = and x = 1 the coexistence region 
is now also rather broader, since for the limiting values of x 
we now recover a system with a Schulz distribution of rod 
lengths, rather than a monodisperse system. 



ist in the phase diagram, the two peaks of the parent 
length distribution must, in some sense, be sufficiently 
well "separated" . This notion is supported by an analy- 
sis of the large q limit in App. There we show that, 
for large q and values of x of order 1/q, there is an N- 
N region and a corresponding three-phase region for all 
as, up to the maximum value as that we consider; cor- 
respondingly, one easily sees that in this limit the length 
distribution always has two separate maxima. 

We now proceed to a closer investigation of the link be- 
tween the presence of three-phase equilibria and bimodal- 
ity of the length distribution, using our second approach 
explained above: we fix q and x so that the three-phase 
coexistence occurs in the bidisperse limit as = 0, and 
then find the value of as at which the three-phase region 
disappears, comparing with the value where the parent 
changes from bimodal to unimodal. 

Figs, [ij and |l^ show for 5 = 12 and two different val- 
ues of X that the three-phase region indeed disappears on 
increasing as above a certain value. This threshold value 
of as is smaller for larger x, showing that in the {po,x) 
phase diagram the top corner of the three-phase region 
shrinks towards smaller values of x as as is increased. 
The x-value marking the bottom end of the three-phase 
region will also depend on as- One would suppose that 
it increases with as until the x-values for the top and 
bottom eventually coincide and the three-phase region 
disappears. However, since in the phase diagrams we 
considered the lower end of the three-phase region is lo- 
cated at extremely small x-values we were not able to 
confirm this explicitly. 
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FIG. 12. Phase behaviour of a bimodal system with q — 12, 
X = 0.015. Shown is the density po of the parent at which 
phase transitions occur, against (on the vertical axis) the 
width as of the two peaks of the bimodal parent. 
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FIG. 13. As Fig. ha but with q = 12, x = 0.03. 



Repeating the above calculations for a range of values 
of X, we arrive at the dashed curve in Fig. |lj: at a given 
X, we have three phase coexistence for values of as below 
the curve. Conversely, horizontal sections through the 
graph give us the extent in x of the three-phase region 
for given as ■ Since the three-phase region cannot extend 
all the way to a; = 0, this shows that the dashed curve 
must eventually bend down towards the horizontal axis as 
X decreases and meet the axis at some small but nonzero 
value of X. We have not been able to detect this part of 
the curve, however, since the x-values concerned are too 
small for our numerics to be accurate. 

We can now analyse whether there is a correlation be- 
tween three-phase coexistence region and bimodality of 
the parent length distribution. For values of x and as 



below the solid curve in Fig. |14[ the parent is bimodal; 
above the curve we have a unimodal length distribution. 
The comparison of the regions under the dashed and 
solid curves suggests that bimodality may be required 
for three-phase coexistence to occur (except possibly in 
the region of very small x, where our data on three-phase 



coexistence are incomplete), but does not necessary en- 
tail it. 

One might therefore look for additional properties of 
the length distribution as predictors of three-phase coex- 
istence, e.g. the overall polydispersity a. Eq. (^0|) showed 
that for a bidisperse system the overall polydispersity a 
vanishes at x = and x = 1 and has a maximum for 
X = l/(q -I- 1); the values of x where three-phase co- 
existence occurs in Fig. ^ are smaller but at least of 
the same order as those where a is maximized. Also, 
the maximum polydispersity in a bidisperse distribution 
with a given length ratio q increases with q, in line with 
the general trend that separation into three phases re- 
quires sufficiently large q. Nevertheless, it is clear that 
the trend of three-phase separation occurring preferably 
for systems with large a cannot be universally true. In 
fact, it is easy to see this from Fig. |lj. When we increase 
the width erg of the two peaks in a bimodal distribu- 
tion the overall polydispersity a also increases; but this 
actually suppresses three-phase coexistence rather than 
enhance it. 

Loosely speaking, our results show that, in order for 
the system to be able to separate into three phases, we 
need a sufficiently large disparity between long and short 
rods in the system, i.e. large q; but at the same time 
the concentration of short and long rods must not be 
too similar. The second condition would account for the 
fact that the three-phase region develops only at small 
values of x, and disappears as as is increased. If this 
interpretation is correct, then we should be able to see 
three-phase coexistence even for a unimodal rod length 
distribution, as long as it contains a sufficiently large 
number of rods that are much longer than the average. 
A good candidate for this is a log-normal distribution, 
which has been shown already in the context of polymer 
solutions to present some interestingJeatures caused by 
the presence of very large particlesEllEj. We will pijesent 
our analysis of this case in a separate publication^. 



V. CONCLUSION 

We have analysed the phase equilibria of the 7^2 On- 
sager model of hard rods with length polydispersity. The 
model is defined by a truncation of the angular depen- 
dence of the excluded volume interaction of the Onsager 
theory after the second Legendre polynomial 7^2 (cos 0). 
Within this model we derived the exact phase equilibrium 
conditions, but these are still rather difficult to solve nu- 
merically. We therefore exploited the fact that the model 
is truncatable - the excess free energy is a function of 
only two moments (pi and P2) of the density distribution 
- and used the moment method, from which numerical 
solutions for phase equilibria can be obtained with well- 
controlled accuracy; the onset of nematic ordering from 
the isotropic side is found exactly by construction of the 
moment free energy. 
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FIG. 14. Shown is, for a mixture of two Schulz distribu- 
tions with length ratio q = 12, the value of as at which we 
lose the bimodality of the parent (solid), together with the 
maximum value of as for which we have a three-phase region 
for the given value of x (dashed). The parent is bimodal in 
the region below the solid curve. Note that, from theoretical 
considerations, the dashed curve has to bend over as x de- 
creases and eventually meet the horizontal axis at some small 
but positive value of x. 



For a fully polydisperse case with a unimodal (Schulz) 
length we found some of the common features of poly- 
disperse systems. In particular, we observed a strong 
fractionation (Fig. ^ and a pronounced broadening of 
the coexistence region (Fig. |l|) with increasing polydis- 
persity; the latter is defined as the normalized standard 
deviation of the length distribution. Fractionation effects 
were strong enough to lead to a crossing of the cloud point 
curves with the corresponding shadow curves in the num- 
ber density representation (Fig. |l|), but not in the volume 
fraction representation (Fig. H). This implies that while 
the nematic phase always has a larger rod volume frac- 
tion than the isotropic, it can actually have a smaller rod 
number density. 

For the unimodal length distribution we did not ob- 
serve any three-phase TN-N or two-phase N-N coexis- 
tence, or re-entrant behaviour of the nematic phase, all of 
which had been found in the Onsager model with bi- and 
tridisperse length distributions. In order to understand 
whether this was due to the introduction of a continuous 
and unimodal length distribution, or to the truncation 
that defined the V2 Onsager model, we next studied the 
bidisperse case within the V2 model. As the length ratio 
q increased, we indeed found (Fig. || and 0) the develop- 
ment of a re-entrance of the nematic phase, and eventu- 
ally also I-N-N and N-N coexistence (Fig. ||) for larger q. 
Thus all the qualitative features of the bidisperse phase 
diagrams of the full Onsager theory are recovered. The 
only exception is the presence in the V2 model of a critical 
point that closes the N-N coexistence region at hidi den- 
sity, while for the full Onsager theory it is knowrO that 
for large q such a critical point does not exist. We ex- 
plained this difference in terms of the different behaviour 



of the excluded volume interaction in the limit of strongly 
ordered rods, and argued that it would persist for any 
truncation of the full Onsager theory involving only a 
finite number of Legendre polynomials. 

Finally, in order to interpolate between the different 
behaviours in the unimodal and bidisperse cases, we in- 
troduced polydispersity into the latter system, by con- 
sidering mixtures of two Schulz distributions of width 
peaked at different rod lengths; the bidisperse system is 
recovered for trg — > 0. For small as, i.e., small "broad- 
ening" of the bidisperse "peaks" , the phase diagrams re- 
mained essentially as in the bidisperse limit, while larger 
values of trg lead to a disappearance of the I-N-N and N- 
N coexistence regions (Fig. |l^). It is natural to suppose 
that this change is related to the fact that the length 
distribution changes from bimodal to unimodal as as is 
increased. However, a closer investigation of this point 
showed that bimodality, although correlated with the oc- 
currence of I-N-N and N-N coexistence, is not sufficient 
to guarantee that such coexistence will occur. Indeed, 
this follows already from the fact that in the bidisperse 
system (which is always bimodal) a minimum value of 
the length ratio q is needed for three-phase coexistence. 
Our results suggest that, also in the polydisperse case, 
a large ratio between the typical (e.g., mean) values of 
the two Schulz length distributions is required. However, 
the result that three-phase coexistence generally only oc- 
curs for small values values of the number fraction x of 
long rods suggests that in addition the length distribu- 
tion must be rather asymmetric, with fewer long rods 
than short ones. When the width as of the Schulz distri- 
butions is increased too much, this (here only vaguely 
defined) "asymmetry" is lost, and the three-phase re- 
gion disappears. In order to assess whether this intuition 
is correct, we are currently investigating the phase be- 
haviour for log-normal length distributionfS; while uni- 
modal, these contain rods of very dissimilar lengths and 
are also strongly asymmetric. The intuition developed 
above would lead one to expect three-phase coexistence. 

We did not investigate in detail the conditions for oc- 
currence of the re-entrance in the nematic phase bound- 
ary. This phenomenon seemed rather robust compared 
to I-N-N and N-N coexistence, although in the unimodal 
(Schulz distribution) case it was again absent. In the 
bidisperse case the re-entrance appeared at lower val- 
ues of q than the three-phase coexistence, and also did 
not require the asymmetry x between long and short rod 
number densities to be quite so extreme, extending e.g. 
almost up to x = 0.3 for g = 12 (Fig. ^). In the case 
of a mixture of Schulz distributions case it also survived 
to larger values of the width as of the two peaks of the 
distribution. An intriguing observation is that, for the 
bidisperse systems, the re-entrance begins (starting from 
small x) very close to the value of x = l/(g-|-l), where the 
polydispersity a of the system is maximal; see Eq. (|4C|). 
This is also the value of x where the volume fractions 
occu pied by the short and long rods become equal; see 
Eq. ( AlOp . We have not been able to find a simple phys- 
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ical explanation for this observation. 

In summary, our study of the 1^2 Onsager model shows 
that the details of the rod length distribution can have 
profound effects on the phase behaviour. In particular, 
while bidisperse systems can show a re-entrant nematic 
phase boundary as well as I-N-N and N-N coexistence, 
these features are absent from the phase diagram for 
a unimodal length distribution; they also disappear as 
one interpolates between the bidisperse and unimodal ex- 
tremes. Put differently, the more "exotic" aspects of the 
phase behaviour of mixtures of hard rods are not ex- 
pected to be observed in systems with only a moderate 
spread of rod lengths around a mean value, but require 
sufficiently wide (and, as we argued above, asymmetric) 
rod length distributions. Our comparison of the predic- 
tions of the V2 model and the full Onsager theory for the 
bidisperse case suggests that these conclusions are not 
artefacts of the truncation used to the construct the V2 
model. Rather, we expect that they would qualitatively 
also be found in an analysis of the effects of length poly- 
dispersity in the framework of the full Onsager model; 
such an analysis remains a challenging problem for fu- 
ture work. 

We comment finally on the presence of the N-N crit- 
ical point within the bidisperse 7^2 Onsager model. As 
discussed above, the theoretical work of van Roij and 
MuldeiEJ shows that for sufhciently large ratios q of rod 
lengths such a critical point will not occur in-tJie full On- 
sager theory; earlier reports to the contraryllj appear to 
be artefacts ot-numerical approximations. However, the 
results of Ref.E3 show only that the N-N coexistence re- 
gion is open, i.e., extends to arbitrarily large densities, for 
q above « 3.17. The possibility remains that for a small 
range of q-values below 3.17 an N-N coexistence region 
would already exist, at moderate densities and closed off 
by a critical point; the topology of the phase diagram for 
such q would be akin to the 1^2 Onsager phase diagrams 
that contain N-N and I-N-N coexistence regions. As q 
increases, the critical point would then have to move to 
larger and larger densities so that the N-N region even- 
tually becomes open at g « 3.17. It would be worth 
revisiting the bidisperse Onsager model to confirm this 
scenario. 

Acknowledgement: PS acknowledges financial support 
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APPENDIX A: THE STRONG NEMATIC 
ORDERING APPROXIMATION 

In this appendix we describe an approximation that is 
useful for analysing the phase diagram of the V2 Onsager 
model at high densities. This strong nematic ordering 
(SNO) approximation becomes exact in the limit ^ 00 
and allows one to show formally that there can be no N-N 
coexistence in this limit. For finite po it is approximate 
but nevertheless useful for understanding some limits of 



the phase behaviour, e.g. for large values of the rod length 
ratio q in bidisperse and bimodal length distributions. 
As discussed in Sec. |VB, it is in the regime of strong 
nematic ordering (small angles 9) that the V2 Onsager 
model and the full Onsager theory differ most strongly, 
so the analysis here focusses on the properties that are 
particular to the V2 Onsager model. 

The SNO approximation is based on the assumption 
that all rod angles in any nematic phases are reasonably 
small, so that we can approximate cos^ ~ 1 — 6^ and 
hence 



7^2 (cos ( 



1 - -i 



(Al) 



From Eq. ([l9|), the orientational distributions are then 
Pi{0) ~ exp[— |^C2P2^^] and so the SNO assumption of 
small angles 9 corresponds to p2 being large enough to 
have lp2 ^ I for all lengt hs I t hat contribute significantly. 
Using the approximation ( |AlD , we can get explicit expres- 
sions for the length dependent chemical potentials. The 
relevant angular integral can be evaluated as 



d9 



JC2P2'P2 



tt/2 



d9 sin 6* €''=^''^"(3/2)''=^''^^" 



exp (/C2P2) 
3/C2/52 



(A2) 



where the symmetry of the integrand under 9 ^ tt — 9 
has been exploited and corrections that are exponentially 
small in lp2 have been neglected. The chemical poten- 
tials (21) thus become 



p.^{l) - lnp{l) + l{cipi ~ C2P2) + H3IC2P2) (A3) 

for a nematic phase, whil e for an isotropic phase {p2 = 0) 
the angular integral (|A2| ) equals unity and the chemical 
potentials are 



pHl)=liip{l) + lcipi 



(A4) 



Using the SNO assumption ( [Al[ ) we can also evaluate the 
integrals appearing in the self-consistency equation (|l^ ) 
for p2', the latter can then be solved explicitly to give, 
after some algebra. 




(A5) 



Po, 
This 



In the limit of large density po, since pi — {l)po 
the square root approaches unity and thus p2 ^ pi 
is as expected, since the weight functions ( p^|jl^ ) for pi 
and p2 only differ by a factor of P2 (cos 9) which becomes 
equal to one for maximal nematic order, i.e. for 9 0. 
The equality p2 = pi implies that p2 becomes large as the 
density increases; for po ^ 00 it diverges and the SNO 
approximation becomes exact. One can now easily show 
that in this large density limit there can be no coexistence 
between two nematic phases. The osmotic pressure (13 
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becomes H 



Coex- 



Po + 5(ci - C2)pI ~ i(ci - C2)pI. 
isting nematics would therefore need to have the same 
pi, hence also the same p2- But then the excess parts 
of the chemical potentials p{l), Eq. (pT|), would also be 
the same in the two phases. Chemical potential equality 
then also forces equality of the density distributions p{l), 
so that the supposed coexisting nematics turn out to be 
the same phase. 



Using the SNO approximation (A3,A5) we can fur- 
thermore get approximate results for the nematic cloud 
curves, i.e. the phase boundaries at high densities be- 
tween the single-phase nematic region and any I-N or N-N 
coexistence regions. The single nematic phase must have 
the density distribution of the parent, p{l) — poP^^^l), 
with pi — po since the parent distribution has unit aver- 
age rod length; its value of p2 is given by (A5). We need 
to find the value of the density po where this phase first 
begins to coexist with a new isotropic or nematic phase as 
Po is lowered. For I -N coexistence, equality of the chem- 
ical potentials (A3,A4) implies that the isotropic phase 
has density distribution 



with 



pHl) = 3lc2P2Poe'^P^°Hl) 



P = Ci(po - P\) - C2P2 



(A6) 



(A7) 



For given /3 and po, p\ is given by p\ = J dllp^{l) and 
the onset of I-N phase coexistence can thus be determined 
simply by solvi ng n umerically for the values of /? and pp 
that fulfill Eq. ( A7 ) and give equal osmotic pressure (^) 
in the two phases. 

For N-N coexistence one similarly finds that the new 
nematic phase has density distribution 



P^ 



with 



/3 = Ci(po - Pi) - C2(P2 - P2) 



(A8) 



(A9) 



and one can again solve numerically, e.g. for /?, po, Pa 
and pf. (The four conditions are then Eq. ([A^), Pq — 
J dl p^{l), Pi = J dl lp^{l), and the osmotic pressure 
equality; p^ is given by Eq. ( [A5| ) in terms of p^ and p^ .) 

We can obtain in this way approximate nematic cloud 
curves for bidisperse and bimodal parent length distri- 
butions; these reproduce qualitatively the features of 
our numerically exact phase diagrams, including the re- 
entrance of the I-N phase boundary. Whether I-N or N-N 
phase separation actually occurs as density is lowered de- 
pends on which of the calculated cloud curves is reached 
first; due to its overestimate of the nematic ordering, it 
turns out the SNO predicts stable N-N demixing only 
above a larger threshold value of the length ratio q than 
in the numerically exact phase diagrams. The overall 
trends are, however, the sa me as in the nu merically ex- 
act phase diagrams of Sec. IVB and IV C. Increasing q 



increases the size of the re-entrant region of the I-N phase 
boundary; the N-N phase boundary also moves to larger 
densities and becomes stable over a wider range. In- 
creasing the peak width erg, on the other hand, produces 
the opposite effects. Interestingly, the SNO predicts that 
the I-N phase boundary remains re-entrant even for the 
largest values of the peak width, 173 = 1- We were not 
able to check this prediction by calculating exact phase 
diagrams for as = 1, but it does suggest that the re- 
entrance is mainly due to the overall polydispersity of 
the parent, rather than the presence of two well sepa- 
rated peaks in the length distribution. 

On the N-N phase boundaries that we calculate within 
the SNO approximation, a critical point occurs, signalled 
by the vanishing of (3, Eq. (A9), and the coincidence of all 
moments po, pi, p2 in the two phases. This critical point 
is not always observable, however, since it can lie in the 
metastable region inside the I-N phase boundary. This 
can happen even if other parts of the N-N phase bound- 
ary are stable. We find that this situation occurs for a 
narrow range of g, for genuinely polydisperse (bimodal 
rather than bidisperse) distributions. This prediction of 
the SNO approximation may appear peculiar but must 
be expected to be a general feature of systems containing 
more than two types of rods. Our (po, x) phase diagrams 
are then planar cuts through a higher-dimensional phase 
diagram; tielines beginning in the cut plane will connect 
to coexisting phases off the plane. For an appropriately 
chosen plane, all tielines originating from a stable N-N 
phase boundary in the plane may then point "to one side" 
of the plane, with none of them having zero length and 
thus giving a critical point. 

Finally, we use the SNO approximation to analyse the 
phase diagrams for bidisperse and bimodal length distri- 
butions in the limit of large rod length ratio g; this limit 
would be difficult to access using direct numerical cal- 
culations of the phase diagrams that do not employ the 
SNO assumption. We also expect the SNO approxima- 
tion to become more accurate for larger q. This is sug- 
gested by the fact that an increase in q generally causes 
the nematic cloud point to shift to higher density; at 
such higher densities the nematic will be more strongly 
ordered, and thus the SNO approximation better justi- 
fied. While such an argument applies directly to the case 
of the boundary between the nematic region and the I-N 
region, the same cannot be said for the N-N phase bound- 
ary. There the hypothesis of strong ordering is applied 
to both nematic phases in the SNO approach, and there 
is no trivial reason why a nematic shadow phase coexist- 
ing with a denser nematic parent is necessarily strongly 
ordered. Numerically, however, we have found this trend 
to be confirmed. 

It turns out that in the limit q — > 00 the equations 
that we need to solve to obtain the nematic cloud curves, 
i.e. the onset of phase separation into I-N or N-N when 
decreasing the density, become dependent only on the 
scaling variable ^ = xq rather than separately on q and 
the fraction of long rods x. This simplification occurs 
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because terms that are exponentially small in q can be 
neglected for q ^ oo. Intuitively, ^ is for large q simply 
the ratio of the volume fractions of long and short rods. 
This follows from the fact that the ratio of number den- 
sities is x/{l — x), so that the volume fractions have the 
ratio 



qx 



hil-x) 1 



(AlO) 



which tends to ^ for x = (,/q ^ 0. The numerical solu- 
tion of the resulting simplified system of equations shows 
that there is always a stable region of N-N demixing for 
q ^ oo; interestingly, this is true even for bimodal dis- 
tributions with the largest pea k wid th cts = 1 that we 
consider. As explained in Sec. IV C, this is presumably 
due to the fact that in the limit q —>■ oo (with x = ^/q 
and ^ held constant) the parent always has two sepa- 
rate peaks, whatever the value of as- The N-N phase 
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FIG. 16. The N-N phase boundary for bidisperse rod 
length distributions with a range of rod length ratios q, as 
derived from the SNO approximation, with the scaling vari- 
able ^ = xqon the vertical axis. The N-N critical point moves 
towards ^ = for q oo, while the N-N boundary as a whole 
moves towards higher densities in the same limit. 
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FIG. 15. Phase diagram as derived from the SNO approx- 
imation, for a bidisperse system with q —> oo; the scaling 
variable on the vertical axis is f = xq. The straight line near 
the bottom left corner of the plot indicates the onset of N-N 
coexistence, which is stable - since it occurs at higher density 
than I-N phase separation - at sufficiently small ^. The same 
is still true if we move to a bimodal distribution, even for the 
maximal peak widths as = 1 that we consider. 



boundary moves towards larger densities as q increases, 
and in fact becomes remarkably simple as q ^ oo: the 
N-N critical point moves towards = 0, and the part of 
the N-N phase boundary at densities below the critical 
point coUapes onto the horizontal axis. A plot of the 
N-N cloud curve for different q with the scaling quantity 
^ — xq on the vertical axis shows this evolution explicitly 
(Fig. ^6|). Note that in the same representation the I-N 
phase boundary would change rather little for the same 
range of values of q, showing that the relative stability of 
the N-N phase separation (compared to I-N) is enhanced 
for large q. This is consistent with the trend o bserv ed in 
the numerically exact phase diagrams on Sec. IVB. 
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